home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_10_04 / 1004030b < prev    next >
Encoding:
Text File  |  1989-12-31  |  1.5 KB  |  52 lines
  1. /* Listing 5 */
  2. typedef struct {
  3.     double X1, X2;
  4. }      ARG_D_2;            /* vector 2 */
  5. ARG_D_2
  6. tan_2(xi1, xi2)
  7.     double xi1, xi2;
  8. {
  9.     double x2, x, n1, x4;
  10.     ARG_D_2 res;
  11. #include "float.h"
  12. #if FLT_ROUNDS != 1
  13.     #error "rounding mode not nearest; adjust code"
  14. #endif
  15. #if FLT_RADIX !=2 && FLT_RADIX != 10
  16.     #error "code not optimum for accuracy in this RADIX"
  17. #endif
  18. #include <errno.h>
  19. #ifndef errno
  20.     extern int errno;
  21. #endif
  22. #include <math.h>
  23. #define M_PI    3.14159265358979323846
  24.     x2 = (n1 = (xi1 > 0 ? 1 / LDBL_EPSILON :
  25.         -1 / LDBL_EPSILON)) + (x = xi1) / M_PI;
  26.     x -= (x2 - n1) * M_PI;
  27.     if (fabs(xi1 / M_PI) >= 1 / LDBL_EPSILON |
  28.     fabs(xi2 / M_PI) >= 1 / LDBL_EPSILON)
  29.     errno = ERANGE;
  30.  /* now in 1st or 4th quadrant */
  31. #define    c0 33281881.3202530279
  32.     n1 = c0 + (x2 = x * x) * (-15666569.8711211851);
  33.     x4 = x2 * x2;
  34.     res.X1 = x * (c0 + x2 * (-4572609.43103684572) + x4 *
  35.        (131095.887915363619 + x2 * (-968.863245687503149 +
  36.            x2))) / (n1 + x4 * (915701.668921990803
  37.                  + x2 * (-13491.7937027796916)
  38.                  + x4 * 44.4083322286368691));
  39.  /* copy 2 */
  40.     x2 = (n1 = (xi2 > 0 ? 1 / LDBL_EPSILON :
  41.         -1 / LDBL_EPSILON)) + (x = xi2) / M_PI;
  42.     x -= (x2 - n1) * M_PI;
  43.     n1 = 915701.668921990803 - (x2 = x * x)
  44.     * 13491.7937027796916;
  45.     x4 = x2 * x2;
  46.     res.X2 = (x * (c0 + x2 * (-4572609.43103684572)) +
  47.        (131095.887915363619 + x2 * (-968.863245687503149 +
  48.       x2)) * x4 * x) / (c0 + x2 * (-15666569.8711211851) +
  49.             x4 * (n1 + x4 * 44.4083322286368691));
  50.     return res;
  51. }
  52.